      subroutine ornl_diir1d(    in_put, incinp, i0_inp, n_inp,
     $			    iirfil, inciir, i0_iir, n_iir,
     $			    output, incout, i0_out, n_out)
      double precision in_put(*), iirfil(*), output(*)
      integer incinp, i0_inp, n_inp
      integer inciir, i0_iir, n_iir
      integer incout, i0_out, n_out
c-----------------------------------------------------------------------------
c
c   dconv1d  performs a 1D recursive convolution in the time domain :
c	O(j) = (1/F(0)) * { I(j) - Sum[ O(i) * F(j-i) ] }, for i=1,...,(j-1)
c
c-----------------------------------------------------------------------------
c
c   PARAMETERS:
c
c	in_put:	Pointer to FIRST sample of sequence "in_put"
c	incinp:	Increment between two successive values of "in_put"
c	i0_inp: Index of the first element of "in_put" 
c	n_inp:	Number of samples of "in_put"
c
c	iirfil:	Pointer to FIRST sample of sequence "iirfil"
c	inciir:	Increment between two successive values of "iirfil"
c	i0_iir: Index of the first element of "iirfil" 
c	i0_iir:	Number of samples of "iirfil"
c
c	output:	Pointer to FIRST sample of sequence "output"
c	incout:	Increment between two successive values of "output"
c	i0_out: Index of the first element of "output" 
c	n_out:	Number of samples of "output"
c
c-----------------------------------------------------------------------------
c
c   PLEASE NOTE:
c    1-	Please note that the array pointers must all point to the first 
c	element of the array "i0_inp", "i0_iir" and "i0_out".
c	 If "in_put" for example is defined as
c	    dimension in_put(-25:45)
c	Then "dconv1d" must be called with the following parameters
c	    call dconv1d( in_put(-25),1,-25,45, ... )
c
c    2-	There isn't any checking done to see if the IIR filter is minimum phase
c	or not. In the case of a non-minimum-phase filter, the computation is
c	numerically unstable, OVERFLOW may results.
c
c-----------------------------------------------------------------------------
c
c   USAGE:
c	This module compute the result of the convolution in the "output" range
c	padding with zeroes when needed.
c
c	With some precautions it is possible that 
c		the Output sequence   OVERWRITE   the Input one.
c	Then (if i0_out == i0_inp), it is necessary that   i0_iir =< 0
c	Said in DSP jargon: "iirfil" must be ANTI-CAUSAL
c
c	In theory, an input sequence of "n_inp" samples starting at time "i0_inp", 
c	filtered by a sequence of "n_fir" samples starting at time "i0_fir",
c	will result in a new signal of infinite length starting at time (i0_inp + i0_fir). 
c	He we just compute here the values that fall in that range and zero the rest.
c
c	 For example when filtering a sequence of N samples, with a
c	filter of m samples. If one wants only to compute the first N resulting
c	samples, the following call can be used:
c	    call diir1d( f, 0, 1, N,  g, 0, 1, m,  h, 0, 1, N)
c
c-----------------------------------------------------------------------------
c
c   This Fortran Subroutine written by
c	Jean-Pierre Panziera
c	Silicon Graphics Inc
c	September 20, 1991
c
c-----------------------------------------------------------------------------

      double precision tmp, zero, one, alpha
      parameter (zero = 0.0, one = 1.0)
      integer i1_inp, i1_iir, i1_out, k0, k1
      integer i, j, kout, kinp, kiir, n0, n1, jmax
c-----------------------------------------------------------------------------
c
c   First Check Parameters
c
c-----------------------------------------------------------------------------
      if( iirfil(1) .eq. zero ) then
	print *, 'In - diir1d - Error - First sample of IIR filter must be != zero'
	stop
      end if

	i1_inp = i0_inp + n_inp - 1
	i1_iir = i0_iir + n_iir - 1
	i1_out = i0_out + n_out - 1
	k0 = i0_inp - i0_iir

	n0 = i0_out
	if( n0 .lt. k0)		n0 = k0

	alpha = one / iirfil(1)
c-----------------------------------------------------------------------------
c
c   Zero the output at the begining of the computated window
c
c-----------------------------------------------------------------------------
	kout = 1
	jmax = min ( i1_out, n0-1)
	do j = i0_out, jmax
	    output(kout) = zero
	    kout = kout + incout
	end do
c-----------------------------------------------------------------------------
c
c	Compute the convolution
c
c-----------------------------------------------------------------------------
	do j = j, i1_out
	    if( j .le. (i1_inp-i0_iir) ) then
		kinp = 1 + (j+i0_iir-i0_inp) * incinp
		tmp = in_put(kinp)
	    else
		tmp = zero
	    endif
	    k0 = max( j-(n_iir-1), n0)
	    k1 = j-1
	    kout = 1+ (k0-i0_out) * incout
	    kiir = 1+ (j-k0)*inciir
	    do i = k0, k1
		tmp = tmp - output(kout) * iirfil(kiir)
		kout = kout + incout
		kiir = kiir - inciir
	    end do
	    kout = 1 + (j-i0_out) * incout
	    output(kout) = alpha * tmp
	end do

      return
      end
